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SO.  abstract 


The  h-version  of  the  finite  element  method  with  piecewise  linear  appro-Kimation  has 
been  applied  to  solve  a  one-dimensional  model  problem  for  the  HelmholU  equation. 
The  main  practical  purpose  of  the  investigation  is  to  lay  theoretical  ground  for  safe 
"rules  of  the  thumb"  how  to  choose  the  meshwidlh  of  the  FE-model  depending  on  the 
wavenumber.  In  this  conte.\t  we  present  new  results  for  stability  and  error  estimation 
of  the  FE-solution.  Following  the  analysis,  numerical  results  are  discussed.  In  a  second 
paper  we  will  study  the  h-p-method  for  Helmholtz  problems. 
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The  h-version  of  the  finite  element  method  with  piecewise  linear  approximation  has 
been  applied  to  solve  a  one-dimensional  model  problem  for  the  Helmholtz  equation. 

The  main  practical  purpose  of  the  investigation  is  to  lay  theoretical  ground  for  safe 
’’rules  of  the  thumb”  how  to  choose  the  meshwidth  of  the  FEl-model  depending  on  the 
wavenumber.  In  this  context  we  present  new  results  for  stability  and  error  estimation 
of  the  FE-solution.  Following  the  analysis,  numerical  results  are  discussed.  In  a  second 
paper  we  will  study  the  h-p-method  for  Helmholtz  problems. 

1  Introduction 

Boundary  _yalue  problems  for  the  Helmholtz  equation  arise  in  a  number  of  physical  ap¬ 
plications  [DL],  in  particular  in  problems  of  wave  scattering  and  fluid-solid-interaction 
[JF].  If  we  analyze  the  scattering  from  an  elastic  body  embedded  in  a  fluid,  analytical 
solutions  can  be  provided  for  regular  shapes  of  the  body  (like,  egs.,  a  sphere  or  a  cylinder 
[JF]).  Numerical  methods  need  to  be  applied  if  the  body  is  of  general  shape.  Here,  the 
physically  proper  and  numerically  effective  modeling  of  large  exterior  domains  is  the  main 
difficulty.  Most  numerical  solutions  have  been  given  starting  from  the  Helmholtz  integral 
equation  applying  boundary  element  methods.  However,  several  difficulties  are  reported 
from  practical  applications  and  finite  element  techniques  are  used  increasingly  not  only 
for  the  solid  but  also  in  the  fluid  domain  (cf.  [H112],  [Bu]).  In  this  context,  the  numerical 
analysis  of  the  finite  element  method  applied  to  Helmholtz  problems  becomes  of  practical 
interest. 

Analytical  results  for  the  finite  element  solution  of  two  point  value  Helmholtz  problems 
in  one  dimension  with  Robin  boundary  conditions  are  contained  in  [AKS]  and  [DSSSj. 
Proofs  of  existence-uniqueness  are  given  for  the  exact  and  the  finite  element  solution 
(h-version)  and  asymptotic  error  estimates  are  proved  under  the  assumption  that  the 
stepwidth  h  is  small  s.t.  the  magnitude  hk"^  (where  k  is  the  wavenumber)  is  small.  This 
assumption  is  obviously  a  severe  restriction  for  practical  applications  if  the  wavenumber  is 
high.  ’’Rules  of  the  thumb”  usid  in  engineering  analysis  of  acoustic  problems  are  given  in 
the  form  hk  =  const  (cf.  [HHl,  p.Tlj).  Some  initial  experience  from  numerical  experiments 
in  fluid-structure-interaction  had,  however,  shown  that  this  rule  failed  in  the  very  case  of 
high  wavenumbers. 
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Hence  the  following  questions  had  been  the  starting  point  for  the  analysis  presented 
in  this  paper: 

•  Are  the  restrictions  imposed  on  h  and  k  in  [AKS]  and  [DSSS]  indeed  necessary  for 
stability  of  the  F&solution  (or  due  to  technicalities  in  the  proof)? 

•  What  is  the  proper  "rule  of  the  thumb”  for  high  wavenumbers? 

•  what  are  the  numerical  and  "economical”  effects  of  the  h-p-version  compared  to  the 
h-version? 

As  in  [DSSS]  and  [AKS],  we  address  these  issues  on  a  one  dimensional  model  problem. 
A  two-point  boundary  value  problem  for  the  Helmholtz  equation  with  Dirichlet  and  Robin 
boundary  conditions  is  considered.  We  start  with  a  recoUection  of  results  for  existence, 
uniqueness  and  stability  of  the  exact  solution  in  the  strong  sense.  We  then  introduce  a 
variational  formulation  for  the  problem,  show  existence-uniqueness  for  the  weak  solution 
and  compute  the  Babuska-Brezzi  stability  constant.  These  results  form  the  prerequisite 
for  the  main  objective,  the  study  of  the  finite  element  solution.  The  analytical  results  of 
this  study  are  contained  in  section  3.  We  first  (3.1.)  formulate  and  prove  a  statement  of 
existence-uniqueness  for  the  FE-solution  using  the  argument  given  by  Douglas  et  al.  in 
[DSSS].  The  proof  is  outlined  in  detail  in  order  to  keep  track  of  all  restrictions  that  have  to 
be  made  for  h  and  k.  The  essence  of  the  argument  is  that  the  FE-solution  is  quasioptimal 
(w.r.  to  k)  provided  the  magnitude  of  hk^  is  sufficiently  small.  However,  quasioptimality 
is  more  than  what  is  needed  in  practical  application  where  (1)  stability  of  the  discrete 
model  and  (2)  error  estimation  at  finite  range  are  the  main  concern. 

Of  these  two  issues,  we  first  address  stability  and  show  that  on  the  finite-dimensional 
level  the  B-B-constant  is  the  same  as  in  the  original  problem  provided  that  the  magnitude 
of  hk  (!)  is  sufficiently  small.  We  then  turn  to  error  estimation  and  show  that  the 
error  is  bounded  if  hk  and  h^k^  are  appropriately  constrained.  Again  this  is  proved  by 
using  an  assumption  on  hk  only.  This  error  estimate  is  the  quantitative  equivalent  to  the 
observation  (cf.  [HHl])  that  in  general  the  error  of  the  finite  element  solution  is  influenced 
by  a  phase  lag  between  the  exact  and  the  finite  element  models. 

In  the  numerical  evaluation  we  present  results  from  various  computational  experiments, 
applying  and  illustrating  the  main  results  of  our  study.  We  show,  in  particular,  that  the 
restriction  of  hk^  is  indeed  necessary  for  quasioptimality  of  the  finite  element  solution. 

2  The  Reduced  Wave  Equation  in  One  Dimension 

In  this  section  we  prove  existence-uniqueness  of  the  solution  to  the  one-dimensional  re¬ 
duced  wave  equation  with  Dirichlet  and  nonreflecting  boundary  conditions.  We  analyze 
the  cases  u  G  H^(Q,  1)  and  u  G  /f*(0, 1)  separately  and  show  that  different  stability  condi¬ 
tions  apply  for  these  two  cases.  The  construction  of  the  Green’s  function  to  the  problem 
is  essential  to  both  proofs. 

2.1  The  boundary  value  problem 

Let  D  =  (0, 1)  and  let  on  il  the  BVP  Lu  =  —f  he  given: 


u"(i) -1- 4:^«(a:)  =  -/(x) 

(2.1) 

II 

o 

(2.2) 

ti'(l)  —  ifcu(l)  =  0 

(2.3) 
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where,  for  simplicity,  /(i)  €  C'(0, 1)  and  k  =  const.,/:  €  IR,/c  >  0. 

Physically,  if  u  is  the  variation  of  pressure  in  an  acoustic  medium  at  a  fixed  time,  eq  (2.1) 
is  the  equation  of  a  plane  wave  with  (nondimensional)  wave  number 


c 


where  u;  is  a  given  frequency,  L  is  the  measure  of  the  domain  and  c  is  the  speed  of  sound 
in  the  acoustic  medium. 

In  X  =  0  a  Dirichlet  boundary  condition  is  given  (prescribed  pressure);  the  mixed  bound¬ 
ary  condition  in  x  =  1  is  the  one-dimensional  Robin  condition. 

Notation:  As  usual,  we  will  denote  by  Z^(0, 1)  :=  .ff°(0,l)  the  space  of  all  square- 
integrable  complex- valued  functions  equipped  with  the  inner  product 

{v,w):=  f  v{x)w(x)dx 
Jo 

and  the  norm  _ 

llu>l|  :=  ^/(^u,u;). 

By  1),  s  =  1  V  2  we  denote  the  Sobolev  space 

=  |u|u  €  A  dui  €  L^,  i  =  1 . .  .s| 

where  dui  are  the  derivatives  of  order  i  in  the  distribution  sense.  The  norm  of  the  space 
H*(0, 1)  is  defined  as 

=  (ii-iiHii-'ii")’. 

Functions  from  /f'  with  Dirichlet  boundary  data  can  be  measured  equivalently  by  the 
/f^-seminorm 

l«li  :=  ll«'l|. 

Uniqueness  of  the  solution  in  1):  The  BVP  (2. 1-2.3)  has  unique  solution  in 
the  space  //^(0, 1). 

Indeed,  suppose  there  exist  two  solutions  M|  and  to  the  BVP  (2.1  -  2.3),  then  u(x)  =  Ui(x)- 
«2(*)  #  0  is  a  solution  of  (2.1  -  2.3)  with  homogeneous  data  f(x)  =  0.  Then  u(x)  is  a 
solution  in  the  classical  sense  and  the  general  solution  of  eqs  (1)  and  (2)  is  u(x)  =  Csin(x). 
Substitution  into  (3)  then  gives 


Ck{cosk  -  fsin/;)  =  0. 


Since 


I  cos  k  —  i  sin  k\  =  1 


we  have  (7  =  0  and  thus 


u{x)  =  0 

which  is  a  contradiction  and  uniqueness  is  proved. 

The  existence  of  the  solution  is  concluded  from  the  following  construction. 
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Inverse  Operator:  The  Green’s  function  of  the  BVP  (2. 1-2. 3)  is: 


G(x,s)  =  - 


j  (  sinArie'*^;  0  <  i  <  s 


sin  kse'^^]  s  <  i  <  1 


The  solution  u(x)  of  (l)-(3)  exists  for  all  A:  >  0  and  can  be  written  as 

ti(i)  =  f  G(x,s)f(s)ds 
Jo 

=  i^sin/fx  J  cos  ksf{s)ds  +  cos  kx  j  s\nksf{s)ds 

-h  tsinA'x  j  sin  ksf{s)ds^  .  (2.5) 

Furthermore,  integrating  in  eq  (2.5)  by  parts  we  see  that  u  €  C^(0, 1). 

Using  the  Green’s  function  we  now  establish  bounds  of  the  exact  solution  and  it’s 
derivatives  by  the  data  /. 

Lemma  1  Lei  u  G  //^(0, 1)  be  the  solution  to  the  BVP  (2. 1-2.3)  for  given  data  f  € 
I2(0,l)=  //■’(O,!).  Then 


l«!l  <  ^ 


lu'll  < 


IK'll  <  (i  +  WII- 

Proof:  Estimates  (2.6)  and  (2.7)  follow  immediately  from  eq  (2.5). 
The  estimate  for  the  second  derivative  is  obtained  from 

=  /\«")'  =  [if-  kW  =  C  2k^  f  fu  + 
Jo  J  Jo  Jo 


<  ll/ll'  +  2t'||/||  ||«||  +  t'll.ll’ 

where  the  Cauchy-Schwarz  inequality  has  been  applied.  With  eq  (2.6)  we  then  get 

ll""ll’  <  ll/ll'  +  2A-^||/||  i||/||  +  A'"||/|p  =  (1  +  k)^Uf  (2.9) 

which  is  the  desired  result. 

Remark  1:  It  can  be  easily  seen  that  all  aforementioned  results  are  valid  also  for  the 
adjoint  problem  (2.1),  (2.2)  and 


li'(l)  -h  ikxt{l)  =  0. 
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2.2  Variational  formulation  and  weak  solution 

A  variational  formulation  of  the  BVP  (2. 1-2.3)  can  be  obtained  formally  by  multiplying 
eq  (1)  with  the  comple.\  conjugate  of  a  suitably  choosen  testing  function  v.  Then,  after 
partial  integration  and  substitution  of  eqs  (2,3)  we  arrive  at  the  variational  problem 

B{n,v)  =  j  (^u'(x)v'(x)  -  k^u(x)v{x)^  dx  -  iku(l)v{l)  =  Fiv)  (2.10) 

where  ^ 

Fiv)=  I  f(x)v(x)dx.  (2.11) 

Jo 

If,  in  general,  there  are  given  Hilbert  spaces  V*  and  and  u  £  (the  trial  space), 
t>  6  (the  test  space)  then  5  is  a  sesquilinear  form 

B  :  V'xV^  -*  €, 


F  is  a.  functional 


F  :  C 


and  a  function  u  €  V*  is  called  a  weak  solution  of  (2.10)  if 


B(u,v)  =  F(v)  (2.12) 

for  all  V  6  V^. 

In  our  case,  the  natural  choices  for  the  trial  and  test  spaces  are 

,/i  ^  j/2  ^  y  ^  //>(o,])Au(0)  =  O}  (2.13) 

For  test  functions  v  €  11^(0, 1),  the  problem. (2. 10)  is  well  defined  if  the  data  /  lies  at  leeist 
in  the  dual  space 


■I 


I/I-,  sup 


1 


Note  that  the  variational  problem  (2.10)  is  equivalent  to  the  BVP  (2. 1-2.3)  in  the  sense 
that  for  sufficiently  smooth  data  any  weak  solution  of  (2.10)  is  a  ’’strong”  solution  of  (2.1- 
2.3). 


Continuity  of  the  form  B:  Applying  Poincare’s  inequality  we  obtain  elementarily 
the  continuity  estimate 

|5(u,t»)|  <  Co(l:)|u|i|u|i 

with  Co  =  1  -f  A:  -I-  k^. 

Existence-uniqueness  of  the  weak  solution:  We  first  show  uniqueness.  Suppose 
again  that  there  e.xist  two  solutions  U\,U2  €  U  to  the  variational  problem  (2.10).  Then 
«  =  Ui  -  U2  0  is  a  homogeneous  solution;  i.e.  u  6  V’  and  eq  (2.10)  holds  with  F{v)  =  0. 
In  particular  we  have  for  v  =  u: 

P(xi,u)  =  J  (^u'(x)u'(x)  -  k^u[x)u{x)^  dx  -  iku{\)u{l)  —  0. 


F.  nUenburg,  I.  Babitika,  Finite  Element  Solution  to  the  Helmholtz  Equation  ...  6 


Since  the  right-hand  side  is  real  this  equation  can  be  true  only  if 

«(1)  =  0. 


Then  it  follows  from  eq  (2.3)  that 

Vr  €  V’  :  /  uv'dx  =  f  uvclx 

Jo  Jo 


and  hence  for  v  =  x: 

0  =  u(l)  -  «(0)  =  /  u'dx  =  f  uxdx. 

Jo  Jo 

Assume  now  /J  ux’^dx  =  0  for  some  natural  n,  then  partial  integration  yields 

0  - - C  ^ - r-  ux^+^dx. 

n+iyo  (n-|-l)(n-|-2)yo 

It  foUows  by  incuction  that 

0  =  /  ui*di  ,s  =  1,3,5, .. . 

Jo 

Since,  as  a  consequence  from  Muntz’s  theorem  [.4,  p.  45],  the  set 

span  {a:*|.s  =  1,3,5,...} 


is  dense  in  Z,^(0, 1)  wo  conclude  that  u  =  0  .  This  is  a  contradiction  to  the  assumption 
and  uniqueness  is  proved. 

For  the  proof  of  e.xistence  we  observe  that  for  the  form  B  a  Gardings  inequality 

7Ze(e(u,«))-hC||u||-'>||u||f  (2.14) 

(where  ||u||i  =  (||w|p  is  the  //’-norm)  holds  for  C  =  C{k)  =  1  -I- 

We  then  have  (sec*,  e.g..  [J,  p.  194])  the  alternative  statement:  either  there  exists  a  non¬ 
trivial  solution  of  tlie  homogeneous  problem  Z/U  =  0  with  Dirichlet  data  0  or  a  solution  of 
Lu  =  f  with  Diriclilet  data  0  exists  for  every  sufficiently  regular  /.  Since  uniqueness  has 
been  proved  existence  follows.  The  proof  is  completed. 


Remark  2:  As  in  the  strong  case  we  remark  that  e.xistence-uniqueness  holds  obviously 
also  for  the  adjoint  form 

5*(u,t')  =  j  (^u'{x)v'{x)  —  k^u{x)fy{x)^  dx  +  iku{l)v{l). 

Remark  3:  .Ml  statements  made  so  far  hold  also  for  the  Dirichlet  condition  u(0)  =  g. 
In  that  case,  the  set  of  admissible  functions  u  is  given  by 

=  |u  €  //'  A  fi  =  t/| , 

this  set  is  related  to  the  trial  space  by  the  bijcctive  map 


ti  =  u  -f  M* 
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where  u  €  V'  and  «*  is  an  arbitrarily  fixed  element  in  U  (see,  e.g.,  [SB,  pp  16/17]  for 
further  detail). 


Stability  in  .^^-norm  and  Babuska-Brezzi-constant:  In  Lemma  1,  stability  of 
the  exact  solution  has  been  proved  for  the  ’’strong”  case  /  €  ff“(0,l); ti  €  H^(0,1). 
However,  recent  computations  [D]  indicated  that  the  stability  estimates  of  Lemma  1  do 
not  hold  for  the  weak  case  /  €  //~*(0, 1);  u  €  //'(0, 1 ).  Indeed  we  will  show  the  following 


Theorem  1  Let  V'  C  /f^O,  1)  be  the  Hilbert  sjxice  defined  in  eq(2.13).  Then  for  the 
sesquilinear  form  B  ;  h  x  V’  —  C  defined  in  eq(2.I0)  the  Babuska- Drezzi  stability  constant 


7  ;=  inf  sup 

v^V 


is  of  order  more  precisely,  there  exist  positive  constants  C\,C‘2  not  depending  on  k  s.t. 

(2.15) 

Proof:  Let  us  first  proof  the  left  inequality  of  (2.15).  We  will  show  that  for  any  given 
u  6  V  there  exists  an  element  v  €  V  s.t. 

mu,v)\>j\W\\\\v%  (2.16) 

Let  ti  €  V  be  given.  Define  v  :=  u  +  z  where  2  is  a  solution  of  the  problem 

Vu;  €  V  :  B{w.z)  =  k‘^(iv,u).  (2.17) 

The  solution  z  exists  and  is  uniquely  defined.  Furthermore,  z  €  H^{0, 1)  and  is  a  solution 
of  the  BVP 


z".^k^z  =  -k\ 

j(0)  =  0 

c'(l)  =  ikzil). 


Then  z  is  the  Green’s  function  transform  of  /(s)  =  k'^u{s). 
The  proof  proceeds  as  follows;  With  r  =  «  +  .r  we  have 


|5(u,  r)|  >  Ref?(u,  r) 

=  Re  {Biu,u)-\- B(u,z)) 

=  Re  (^B{u,u)-\- B{u,z)  +  k'^(u,u)-  k'^(u,u)^ 

=  ReB{u,u)-i-k^uf  =  \\u'\\\ 


Then,  if  we  show  that 

IIk'II  >  f 

we  have  proved  ineq.  (2.16)  and  the  inf-sup-condition  follows. 

To  obtain  ineq.  (2.18)  we  integrate  by  parts  the  function 

^(a;)  =  k^  f  G{x,s)u(s)ds 
Jo 

=  k^  (^11  {x,  l)rt{\)  -  j  H[.r,s)u'(s)ds 


(2.18) 
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where 


H(X,s):=  r  G{x,t)dt 
Jo 


and  the  integral  on  the  r.h.s.  is  well  defined  since  u  €  //'(0, 1). 
We  now  differentiate 


z'{x)  =  k^  ^Hr(x,l)u(l)-  HAx,s)u\s)ds^ 


where  the  function  //,  is  readily  computed  as 


Hr(x,s)  =  -- 


;^•Ie***  0  <  X  <  s 

iA:se‘**  s  <  X  <  1 


Then 


|x'(x)l  <  k^(\H,ix,l)\\u(l)\  + j\HAx,s)u'{s)\ds^ 

<  A:2(|^.(x,  1)1  + 11^,11)  ||ti'||<2/:||u'|| 

since  obviously  \Hr(x,  1)|  <  |,||/fx||  <  p 
Hence 

lMl  =  ||«'  +  ^11<M  +  ||/||<(l  +  2/:)||u'|| 

or  ^ 

11“'"  ^  ^  t"“'" 

for  k  >  I  and  the  first  part  of  the  proof  is  completed. 

To  proof  the  second  inequality  it  is  sufficient  to  find  some  function  Zo{x)  €  V  for  which 

^  ^  C,  , 

TjT-i'*’'" 


Consider  the  function 


sin  kx 


where  ip  €  C°°(0, 1)  does  not  depend  on  k  and  is  choosen  s.t. 

z,(0)  =  x„(i)  =  x.;(o)  =  x;(i)  =  o 
and  that  for  some  a  >  0,  not  depending  on  k  , 

Noll  >  ct 

(take,  e.g.  (^(x)  =  x(x  -  1)^). 

■Then 

Noll  Q 

and  with  eqs  (2.19)  we  obtain  by  partial  integration 

Vr€  V'  :  P{z,.v)=-  l\z':  +  kh^ 

Jo 


(2.19) 
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Direct  computation  shows  that 

i'o  +  =  s?"  +  2if'(x)coskx. 

Define  ^ 

«(*)  :=  (<(^)  +  k^2o{s))  ds, 

then 

|5(2o,v)|=  u(1)u(1)-  /  u(i)i;'(a-)dx  <  (|u(l)|  +  l|t/||)|u|i. 

Jo 

On  the  other  hand,  from  the  definition  of  u  we  obtain  integrating  by  parts 
|u(i)|  =  jy  ^(^"(5)  +  2(,?'cosfcs^  dsj 

I...  ,,  ^smkx  ,  sinfc5'\  ^  | 

Hence 

ll«ll  <  J  (llv"IU  +  2^'IU 
so  there  exists  a  constant  C  s.t. 

(l«(i)l  +  IHI)<f. 

Consequently, 

Vi-€F:  |5(r„,t;)|<  ||v|, 

and  the  proof  is  completed. 

From  general  theory  p  i  12]  now  follows 

Corollary  1  Let  11  €  //'((I.  1 )  6c  0  so/m/?oji  of  the  variational  problem  (2.10)  with  given 
data  f.  Then  the  .stability  c.>ilimalc 

tw|i  < 


holds  where  C  is  a  generic  constant  not  depending  on  k. 
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3  Finite  Element  Analysis  of  the  Wave  Equation 

Following  preliminary  riefinilions  we  state  appro.ximability  of  the  e.\act  solution  as  a  direct 
conclusion  from  the  approximation  properties  of  the  finite  element  space  and  stability  of 
the  exact  solution  (3.1).  We  apply  an  asymptotic  approach  (|,4KS],  [DSSS])  to  obtain 
statements  of  stability  and  quasioptimal  error  estimates  for  the  model  problem  (3.2). 

In  the  second  part  of  our  analysis  we  study  the  properties  of  the  finite  element  solution  in 
the  preasymptotic  range.  In  a  preliminary  subsection  (3.3)  we  analyze  the  discrete  model 
on  uniform  mesh  and  construct  a  discrete  Green’s  function.  We  then  show  the  inf-sup 
condition  (3.4)  and  state  a  prea.symptotic  error  estimate  (3.5).  The  section  is  concluded 
with  some  comments  (3.6). 

3.1  Approximability  of  the  exact  solution  in  the  finite  element  space 

Preliminaries:  To  solve  the  problem  numerically  with  a  finite  element  method,  the 
inter'''al  Ti  =  [0, 1]  is  divided  into  n  finite  elements 

n 

^  ’  0  =  <  2:1  <  . . .  <  x„  =  1.  (3.1) 

j=i 

By  (3.1)  a  discrete  subset 


•V/,  =  {.rj,j  =  0. 1 . 11}  C  [0.1]  (3.2) 

{finite  element  mesh)  is  given  on  fl  .  Any  function  defined  on  Xh  is  called  a  mesh  function 
and  will  be  referred  to  by  subscript  h. 

We  will  denote  l)y  hj  =  (Xj  )  the  .srcf  of  the  finite  element  #  j  and  define  the  stepsize 

h  of  the  mesh  Xh.  by 

h  =  muxhj.  (3.3) 

j 

In  the  following  we  will  seek  approximate  solutions  of  the  variational  problem  (2.10)  using 
the  Galerkin  finite  element  method  with  piecewise  linear  test  functions. 

More  precisely,  we  define  the  subspace 

5a(0,1)C  //'(0,1) 

as  the  set  of  all  functions  u  €  H^{0,1)  such  that  the  restriction  of  u  to  any  element 
is  a  linear  function.  On  each  finite  element  [xj_i,x^],  a  function  v  6  5/,(0, 1)  is 
written  by  means  of  the  nodal  shape  functions  (SB,  p.38]  A'j^\x),  as 

v(x)  =  i>_,_,.'V|-'’(.r)  -f  VjNj-'\x) 

where  Uj_i.  I’j  a'v  the  nodal  vab'es  of  v  in  x,^\.Xj,  respectively. 

An  admissil)le  function  vj.{.v.h)  €  H'iil)  will  ire  called  a  finite  element  solution  to 
the  variational  problem  (2.10)  if 

1.  u/e  €  .‘>'/.(0,  1)  ; 

2.  7i/e(0)  =  g-. 
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3.  u/e  is  a  solution  to  the  variational  problem  (2.10)  for  all  test  functions  v  € 
where 

Vh  :=  Sh[Q,  1{=  {v  €  Sh(0,  1 )  A  t;{0)  =  0}  . 


Remark  4:  As  in  the  continuous  (cf.  Remark  2)  case,  test  and  trial  space  are  identical 
in  the  discrete  case: 

=  =  Sh[0,l[. 

Approximation  properties  of  5/i[0, 1[:  It  is  a  well  known  fact  that,  for  any  function 
u  6  in  one  dimension  the  best  approximation  on  given  mesh  Xh  is  obtained  by 

the  linear  interpolant  u/  of  u.  Furthermore,  if  «  €  //^(fl),  the  following  statements  hold: 

Lemma  2  Let  u  €  //^(0, 1)  and  ui  £  .Sa(0,  1)  be  the  piecewise  linear  interpolant  of  u. 
Then 

inf  ||w-r||  <  ||m-«/|1  <  ||u"||  (3.4) 

\ir  j 

inf  ||u"||  (3.5) 

veSh  \7r/ 

||«  -  «/||  <  ^^^  l«  -  v|i  (3.6) 


Proof:  see,  e.g.,  (SF,  p.  45]. 

A  statement  of  approximability  is  now  immediately  obtained: 

Theorem  2  Letu  £  //^(fi)  be  the  solution  of  the  variational  problem  (2.10)  -  or,  equaiva- 
lently,  of  the  BVP  (2. 1-2.3)  -  for  given  data  f  £  L^{Q). 

Then  the  for  the  error  of  the  best  appmximalion  in  H  ’  -seminorm  there  holds 

\n-n,U<-{\-¥k)\\f\\. 

n 

Proof:  The  statement  follows  directly  from  Lemmas  1  and  2. 

Remark  5:  The  practical  conclusion  from  this  theorem  is  that,  for  any  given  /,  we  can 
control  the  approximation  error  by  bounding  appropriately  the  magnitude  of  hk.  More 
precisely,  for  any  given  data  /  £  Z.^(0, 1)  and  error  bound  s  >  0  there  exists  ^  >  0  s.t.  for 
hk<6 

inf  («  — uIj  =  <  t. 

Since  the  wave  number  k  is  related  to  the  (noudimensional)  wavelength  A  by 


27r 


the  product  kh  is  a  measure  of  the  number  of  elements  per  wavelength. Hence  the  "rules 
of  the  thumb”  recommending  a  certain  number  of  mesh  points  per  wave  length  do  apply 
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effectively  for  approximability  (of  the  exact  solution)  by  piecewise  linears. 

The  essential  question  is,  however,  if  this  rule  does  also  apply  for  the  finite  element  so¬ 
lution.  Since  the  approximation  property  has  been  established,  the  answer  to  this  question 
would  now  be  given  by  an  stability  estimate  of  the  form 

!t<  -  «/e|i  <  Ct  inf  |u  -  r|i 

where  the  stability  constant  C,,  in  general,  depends  on  k.  Before  turning  to  this  estimate 
we  verify  that  the  variational  problem  (2.10)  is  well-defined  also  an  the  finite  level. 


3.2  Asymptotic  stability  and  quasioptimal  error  estimates 

A  proof  of  existpnfo-iiniqiionc.ss  in  terms  of  asymptotic  error  estimates  of  the  FE-solution 
for  the  one-dimensional  Helmholtz  equation  with  non-reflecting  boundary  conditions  has 
recently  been  given  by  Douglas  et  al.  [DSSS].  Since  we  are  interested  in  the  dependence  of 
error  estimates  on  h  and  k,  specifically  in  the  case  of  large  k,  we  outline  here  the  argument 
from  [DSSS]  (slightly  modified  to  account  for  the  Dirichlet  condition  at  i  =  0)  in  detail 
to  keep  close  track  of  the  constants  involved  in  the  estimates. 

Theorem  3  Let  u  6  //^(0,1)  be  tl  exact  and  ujf  €  5/,(0, 1)  be  the  finite  element  solu¬ 
tions  of  the  BVP  (2. 1-2.3),  resjyectivcly. 

The  finite  element  solution  is  then  uniquely  delermind  by  any  data  f  €  £^(0,1)  and, 
furthermore,  the  following  error  estimates  hold 

||«-«/ell  <  CiC2{\  +  k)-^h^\\f\\  (3.7) 

<  C2{l  +  k)h\\f\\  (3.8) 

with 

r  2 

(1 -2(1-1- A-)^)7r 

and 

C  2  •“  1 

TT  [l-GCfkViHi +kyy 

provided  that  the  siepwidtli  h  and  the  wavenumber  k  aie  such  that  the  denominators  of  the 
consUinis  are  iiositive. 

Proof:  Denote  e  :=  u  -  ujg.  Then  e  lies  in  the  Hilbert  space  V  C  £f^(0, 1)  and, 
consequently  (cf.  remark  3),  there  exists  ^  6  V  s.t. 

Vi'  €  V'  :  B(v,z)  =  (u, e). 


In  particular,  f?(e,.r)  =  (e,e)  for  u  =  e. 

Further  the  error  is  f?-orthogonal  to  the  discrete  test  space  Vh  :=  5/,[0, 1[: 


Vtn  G  Vh  :  B(e,  w)  =  0. 
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Then,  for  all  w  €  Vh, 

lkll^  =  (e,f)  =  B{e,z-zi') 

=  J  e'iz  -  u>y  j f  (-f  -  w)  -  ifce(l)(^(l)-  u>(l)) 

<  II  (-’  -  u>ni  lie'll  +  k^\\z  -  t.11  ||e||  +  k\z{l)  -  n;(l)||e(l)l 

Apply  the  inequality  |v(l)|  <  v/2||t;||i||t>'||»  which  is  true  for  aU  i;  6  to  obtain 

*1^(1) -tn(l)||e(  1)1  <  2^'||  (i -u))'||a  lie'll  j||z-u)||j|le||2 

<  A-^lk  -  HI  lle||  + II  (^-uO'll  lie'll  (3.9) 

where  the  inequality  2a6  <  +  6^  has  been  applied. 

This  gives,  for  all  w  €  V),, 

IMP<2(ll(-'-»-)'lllk'll  +  t'"ll*--HI  W)- 

In  particular  we  may  apply  Lemmas  1  and  2  for  tu  =  r/  G  14  (the  piecewise  linear 
interpolant  of  z)  to  obtain 

M'  <  (iK-- ll<-'ll  +  ‘'’ll---=(IIIWl) 

<  2  (( 1  +  !•) file'll  ll<’ll  +  +  <■•)  Hell’)  ■ 

Divide  both  sides  of  the  inequality  above  by  the  common  factor  ||e||,  then 

||e||<C,(l  +  /:)/i  lie'll  (3.10) 

holds  with 

r  2 

’■  (l-2(H-fc)^)7r' 

under  the  assumption  that  k,  h  are  such  that  the  denominator  of  Ci  is  positive. 

Next,  from  £?-orthogonality  of  the  error  to  elements  from  14  we  have 

f?(  r.  r )  =  P{  e.  u  -  Vf()  =  P( c,v) 


and  hence 

Vr  Gift:  P{e.  r)  =  P(e.  u  -  e). 

Thus,  for  all  i»  G  1  a 

J  e'ef  ~  k^  J  ce  -  ik\c{])\^  =  J  c'{u  -  e)'  -  k^  J  e{u  -  v)  -  74:e{l)(u(l)  -  u(l)) 
and  therefore 

||e'||2  <  A:'^||el|2  +  A'|e(l)|2  +  |(r'||  1|(«  -  r)'l|  +  A’2||e||  |(u  -  r||  +  A|e(l)l  |u(l)  -  u(l)| 
<  A^llcll^  +  2A-||e'l|  ||e||  +  2||c'|l  ||(«  -  r)'||  +  2k^\\e\\  ||«  -  t;|| 
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where  the  terms  in  i  =  1  have  been  estimated  as  in  (3.9).  We  now  use  the  £-inequality  to 
f^et  the  estimates 


2/^-lk'IIIHl  < 

2|le'||||(u-i>)'|l  < 

2A:2  1H|||u-u||  < 


L\le'f  +  4k^\\ef 

llle'lP  +  4|Kw-i.)T 

k^  Help  + -  ull^ 


Introducing  these  estimates  into  the  inequality  leads  to 

Vt;€l4:  IkT  <  6k^ef+\\\er  +  m^-v)T  +  k>-v\\^- 


(3.11) 


Then,  using  the  intermediary  result  (3.10)  and  the  appro.ximation  results  from  Lemma  2 
for  V  =  u/,  we  get 


6kHl  +  ^  ll/ll'- 


and  hence 

+  <  (I)  (i+ (s)')*''(i  +  *)B/ll 

and  the  statement  of  the  theorem  follows  .  The  proof  is  completed. 


Remark  6:  For  the  denominator  of  C'2  to  be  positive,  the  magnitudes  of  {hky,  h^k^ 
and  h^k‘*  need  to  be  small.  The  term  (hk/2it)'^  in  the  numerator  can  then  be  omitted. 


Let  us  state  as  a  corollary; 

Corollary  2  With  the  assumptions  of  the  theorem,  the  estimate 


holds  for 


rvith 


l^ex  “  «/e|i  <  Cs  inf  |u  ex  v|i 

v€»h 


C,:= 


C,  := 


-  dC^k^h^l  +  k)^y 
2 


(1-2(1  +  A’)^^);r 
Proof:  Introduce  cq  (3.6)  from  Lemma  2  to  (3.11). 


(3.12) 


Remark.  7:  Note  that,  except  for  the  final  estimates  in  terms  of  ||/||,  the  proof  of 
the  theorem  is  valid  also  with  the  weaker  assumptions  u  G  /f’(0, 1),/  €  /f~^(0,l).  In 
particular  we  obtain  the  statement  of  quasioptimality  from  the  corollary.  However,  the 
assumption  that  k^h  be  small  is  essential. 


F.  Ililenburg,  /.  Dabuska,  Finite  Element  Solution  to  the  Helmholtz  Equation  ...  15 


Hence  if  h  and  k  fulfill  the  a.ssumptions  of  the  theorem,  the  finite  element  solution 
behaves  effectively  like  the  be.sl  appro.ximation  (i.e.  C\  can  be  replaced  by  some  ab.colute 
constant  not  depending,  on  A  )  and  the  "rules  of  the  thumb"  apply  for  the  FE-solution. 

The  assumptions  of  the  theorem  imply,  however,  that  the  magnitudes  of  hk,  hk^  and 
hk^  have  to  bo  bounded  by  sufficiently  small  magnitudes  (cf.  also  assumptions  in  [DSSS, 
p.  177]).  The  theorem  (and  the  corallary  as  well)  then  states  that,  with  these  restrictions 
to  h  and  k,  the  error  of  the  fe-so!ution  is  quasioptimal.  While  this  is  the  desired  result  it 
is  achieved  at  high  cost  if  k  is  large  and  the  stepsize  h  must  be  choosen  s.t.  the  magnitude 
hk^  is  small. 

At  this  state  of  our  investigation,  it  is  not  clear  whether  the  assumptions  of  the  theorem 
are  due  to  technicalities  of  the  proof  or  really  inherent  to  the  problem  considered. 

The  second  and  more  important  question  is  whether  the  assumptions  of  the  theorem 
are  indeed  necessary  to  bound  the  discretization  error  by  some  finite  magnitude  (like,  egs., 
a  given  tolerance  for  the  relative  error).  The  following  simple  computation  indicates  that 
this  is  not  the  case  for  high  k.  Indeed,  let 

hk^  <  ft 


for  some  ft  >  0.  Then  h  <  a/k^  and 

|l'  -  M/ell  <  C’zd  +  /v)pll/ll 

hence  the  error  estimates  of  the  theorem  tend  loivards  0  (while  they  have  only  to  be 
bounded  for  practical  purposes)  as  k  is  increa.sed. 

We  will  state  stability  under  weaker  assumptions  and  give  more  appropriate  error  esti¬ 
mates  in  a  preasymptotic  analysis  using  a  discrete  Green’s  function  approach  on  uniform 
mesh. 


3.3  Preasymptotic  analysis:  Preliminaries 

Global  FE-equations  and  discrete  fundamental  system;  Let  in  the  following  the 
FE-mesh  be  uniform  with  /j  =  i.  After  assembling  the  local  equations  (2.10)  and  mul¬ 
tiplying  the  whole  set  by  h.  we  arrive  at  a  set  of  linear  equations  for  the  mesh-function 
«/,  =  Ufelx^: 

LhUh  =  Th  (3.13) 

where  the  discrete  operator  L/,  can  be  written  as  a  n  x  n-tri diagonal  matrix 


with 


and 


■  2,S'(/) 

/?(/) 

/?(/) 

2.S(/) 

m 

11(1) 

2S(t) 

mt) 

- 

n(t) 

S(()- 

/?(0  =  - 

'-1 

.  .S’(/) 

/  =  hk. 


(3.14) 
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The  right  hand  side  r/,  is  a  mesh  function  obtained  from 


Vj: 


fix)  N2^(x)dx  + 


fix) 


(3.15) 


Remark  8:  As  noted  before,  the  product  i  =  kh  is  a  measure  of  the  number  of  elements 
per  wavelength  (of  the  exact  solution).  In  particular,  if  the  stepwidth  is  such  that  t  =  f 
for  integer  /  then  exactly  I  elements  are  placed  one  one  half-wave  of  the  exact  solution. 


For  later  use  we  introduce  difference  notation  as  follows:  Given  a  mesh  function  u  =  m 
defined  on  A'/,  we  will  denote  left  and  right  differences,  respectively,  by 

ji  u(ar,)  -  «(ij_i)  u(xi+i)  -  u(x,) 

a  u  := - 7 - ,  u  u  .=■ - 7 - . 

hi  ht+i 

In  the  linear  space  of  mesh  functions,  an  inner  product  in  /.^-analogy  is  defined  by 

n 

ifh.ah)h  =  *  H  fjS]- 

j=i 


We  will  denote  the  discrete  /!,^-norm  defined  by  this  inner  product  by  ||  •  ||.  The  discrete 
analogon  to  the  //’-seminorm  is  given  by 

1  =  1 

Note  that  for  any  piecewise  linear  function  u  with  nodal  points  on  Xh. 

|«li  =  lItt'lP  = 


i.e.  the  discrete  and  exact  /f '-norms  are  identical.  We  will  use  the  discrete  Dirac  symbol 


defined  as 


1  if  » =  i 

- 

0  if  i^  j 


Discrete  wavenumber  and  Green’s  function:  The  fundamental  system  of  eq 
(3.13)  is 

Fn  =  €  {j/mj  =  0, 1 . n}}  (3.16) 

where  k'  is  a  parameter  to  be  yet  determined. 

To  this  end,  we  solve  any  of  the  ’’interior"  equations  in  the  point  Xj  =  j/n,  1  <  j  <  n: 

^  25(/)e‘‘''-''*  +  fl(/)c*'-'<^+'>'*  =  0.  (3.17) 


A  = 


With 
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eq  (3.17)  has  the  solutions 


'SHD 


ii(i)  V 


-1  = 


(♦)  complex  conjugate  if 


(♦♦)  real 


if 


1^ 

1S« 


<  1 


(3.18) 


(<]|  -  ^ 


From  the  definition  of  A  we  see  that  the  discrete  wave  number  k'  is  either  real  (in  case  (*)) 
or  pure  complex  (case  (**)).  Physically,  ca.se  {*)  describes  a  propagating  wave  whereas 
case  (**)  describes  a  decaying  wave  [Hill].  For  sufficiently  small  h  (more  precisely,  for 
h  <  y/T2/k)  one  obtains  always  the  complex  conjugate  solution  of  case  (*). 

The  discrete  wavenumber  k'  can  be  formally  computed  in  terms  of  t:  From  eq  (3.18), 
case  (♦),  we  get 


and  hence 


//  1 

k  =  —  arccos 
k 


Consider  the  Taylor  expansion 

k'h  =  arccos 


21 


Hence,  for  fixed  k. 


A  '  =  A-  - 


k^h^ 

24 


S(t) 

m‘ 

(3.19) 

(  sm 
\  R{i)J' 

,  V 

(3.20) 

D 

► 

-  -I-O(A-V) 

(3.21) 

Once  the  discrete  wavenumber  has  been  computed,  a  discrete  Green’s  function 
Gh{x,s)-,  X  =  Xh,^  =  Sh  can  be  constructed.  We  give  next  a  brief  outline  of  this  construc¬ 
tion  referring  to  [Sa]  for  details. 

Similarly  to  the  continuous  case,  we  require  that  the  r.h.s.  of  the  linear  system  (3.13) 
is  mappped  to  the  discrete  solution  of  this  system  by 

Uh{x)  =  (6'/,(a-,s).r/,(s))/,. 

We  accordingly  seek  the  discrete  Green’s  function  in  the  form 

Cl  sin  k'x  X  <  s 


Gh(x>s)  =  < 


(3.22) 


(  C2(  sin  k'x  -f-  cos  k'x  ) 


s<x<l 


where  C],  C2  are  functions  of  s  and  the  constant  .4  is  determined  from  the  discrete 
equation  in  the  nodal  point  .r„  =  1  as 

^  _  sin  k'  cos  A'(  l}( i sin^  k'h  -  t^)  —  itR{t)  sin  k'h 
/?( i  cos^  k'  sin^  k'h  -|-  t  ^  sin  '^  k' 

Since  o(.r)  :=  .sin  k'x  and  ii(x)  :=  .4sinAr'.T  -f  cos  k'x  are  fundamental  solutions  of  eq. 
(3.13).  we  can  prove  by  discrete  Green’s  formula  ([Sa,  pp.  120/121])  that 


A(.Tj)  =  (d-'a)S  -  cx(d^j3)  =  const  = 


sin  k'h 
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Using  this  formula  we  find  from  the  condition 


=  -f 

the  unknown  coefficients  Ci ,  Cj  to  be 

Pis) 


Ciis)  = 
C2i^-)  = 


Ms) 

q(^) 

A{s) 


Introducing  these  results  into  eq  (3.22)  the  discrete  Green’s  function  is 

sin  k'x{A  sin  k's  +  cos  k's)  x  <  s 
sin  sin  k'x  +  cos  k'x)  s  <  x  <  I 


Ghix,s)  = 


1 


h  sin  k'li 

and  the  discrete  solution  Uh(xh)  =  Gh{xk-,Sj)rh{sj)  becomes 


(3.23) 


h  sin  k'h 
for  0  <  /  <  n. 


if'  ”  ” 

— fth  I  ^  XI  ^  XZ 

j=i 


j=/+i 


j=i 


(3.24) 


Remark  9:  A  straightforward  asymptotic  analysis  of  the  discrete  solution  shows  that, 
for  h  — ►  0  the  coefficient  A  converges  to  i  and  «/,(x/,)  converges  to  the  exact  solution  u(i) 
as  given  in  the  previous  section. 

Remark  10:  Using  eq  (3.19)  the  constant  .4  in  the  Green’s  functions  (eq  3.22)  can  be 
simplified  to 


A  = 


sin  k'  cos  k'  +  iy/T2\/l2  — 


12-/^  cos‘  k' 

Obviously  |/1|  is  bounded  independently  of  k  for  t  =  hk  <  a  <  v/l2. 


(3.25) 


3.4  The  inf-sup-Stability  Condition  for  the  Finite  Element  Solution 

In  this  subsection,  we  will  compute  the  IJabuska-Brezzi  stability  constant  of  finite  element 
solutions  on  uniform  mesh  using  the  discrete  Green’s  function.  E.xistence-uniqueness  of 
the  FE-solution  then  follows  under  weaker  (compared  to  the  proof  outlined  in  the  previous 
subsection)  assumptions  on  h  and  k. 

Stability  of  the  finite  element  solution  and  discrete  B-B-constant:  The  sta¬ 
bility  investigation  of  the  form  B  on  the  finite  level  is  proceeded  in  close  analogy  to  the 
infinite-dimensional  case  as  considered  in  section  2.  Namely,  we  will  prove 
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Theorem  4  Let  14  :=  5’/,[0, 1[€  1)  and  B  :  14  x  V'/,  —  C  be  the  sesquilinear  form 

defined  by  equation  (2.10). 

Then,  if  the  stepwidth  h  <  ^  (or,  resjHctively,  hk  <  1),  the  Babuska- Brezzi  stability 
condition 


inf  sup 


|g(u.i>)| 


=  7A  >  0 


(3.26) 


holds  and  there  exist  positive  constants  C\  and  Cz,  not  depending  on  k  or  h  s.t. 


£i 

k  ■ 


Proof:  The  line  of  proof  is  similar  to  the  infinite-dimensional  case:  we  will  show  that  for 
any  given  «  €  14  there  exists  some  n  €  14  s.t. 


|B(f.»)l>  jMl  Ill'll. 

Let  hence  u  £  14  be  given  and  define  v  :=  11  + z  where  z  £  1 /,  is  a  solution  of  the  variational 
problem 


V«-€14:  B(u\z)  =  k^(w,u).  (3.27) 

Since  14  is  a  Hilbert  space,  the  .solution  of  (3.27)  exists  and  is  uniquely  defined. 

As  in  the  continuous  ca.se,  we  will  now  prove  that 

l-li  >  f  Ivli 

using  the  Green’s  function  representation  of  z: 

n 

=  ShiXi)  =  h  ^  Gijrj  (3.28) 

i=i 

where 

Gij  :=  Gh(xi,Sj);  rj  :=  r^isj). 

Summation  by  parts  in  eq.  (3.28)  yields 

n 

'i  =  ffi„r„  -  HiiTf,  -  h^  Hijdh  (3.29) 

j=i 

with 

£>•'  n,.  =  r/y.  y  =  1 , . . .  R  -  1.  (3.30) 

Since  the  mesh  function  Jl  is  defined  by  eq  (3.30)  up  to  a  constant  we  are  free  to  choose 


//„  =  0. 

Let  us  now  take  the  left  differences  of  Zh  in  some  fixed  point  i  =  /: 

d‘z  =  d‘lI.„r„-hj2<i‘H.jd^r. 

j=i 


(3.31) 
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Then,  applying  the  Schwarz  inequality,  we  obtain  the  estimate 

\d‘z\  <  M'//.„||r„|  +  ||//x|||r|, 

<  (|rf'//.n|  +  ||//r||)|rl,. 

The  right  hand  side  of  the  variational  problem  is  by  direct  computation 
fj  =  ^k^h^(uj.i  +  4uj  +  Uj+i ) ,  j  =  1 . . .  n  -  1 

hence 

|r|i  <  Ch'^k^\u\i 

where  C  is  a  constant  of  order  1. 

We  now  turn  to  estimation  of  the  magnitude  \d‘H.„\  +  ||//r||. 

Prom  eq  (3.30)  we  obtain  after  summation  over  j: 

j-\  j-J 

//,,  -  Hi,  =h^  D‘H,.  =  h^  Gu 

/=!  /=1 

and  consequently,  since  Hu  =  0. 

j-i 

H.j^hJ^Gu  . 

/=i 

Taking  left  difTcrences  wo  obtain 

(PH.j  =  h'^d'G.i 
/=! 

The  derivatives  (as  loft  differences)  of  the  discrete  Greens  function  are 


d'  sin  k'xh  ( A  sin  k'si  +  cos  k'si)  Xh  <  s/ 

< 

sin  k'si  (Ad'  sin  k'xh  +  d'  cos  k'xh)  Xh  >  s/ 


We  substitudo 


„  .  2  fk'li  .  \  .  k'h 

d  sill  A  .I'l,  =  -  ros  I  — (2/  -  1)  sin  — — 

//  V  2  /  2 

d‘  fos A'.t/,  =  —J  sin  ~  ^0 


to  obtain 


(PG.i  = 


/j2  cos  ^ 


cos  ^^(2*  -  1))  (/I  sin  k'si  +  cos  k'si)  i  <  I 

sin  k'si  ^/1  cos  ^^(2*  -  1)^  -  sin  (^(2j  -  1)^^  i  >  I 


(3.32) 


(3.33) 


(3.34) 


(3.35) 


(3.36) 


(3.37) 


F.  Ihlenburg,  1.  Babuska,  Finite  Element  Solution  to  the  Helmholtz  Equation  ...  21 


Then 

Isl 


I  = 


^cos(^(2.-l)j  (^A'^m.k'hl  +  ^cosk'hl^ 

^sinib'/t/  ^i4cos  “  0^  -  sin  ~  ^0)) 

t  A'/,,.,.  .A  Z'.  .  (iVj)k'h  .  (,Vj  +  l)i'A 

V  I”  *  /  V  - 2 - ®‘" - 2 - 


1 

h- cos  —-sin  — 


+  l)k'h  .  ik‘h  .  {i  +  \)k'h'\ 

— - sin  — —  sin - - -  I 

2  2  2) 


.  (^,os(^(2,-l))-sin(^(2.-l)))) 


r>i 


-  h-s\nk'h 

since  |/1|  and  hence  the  e.vpression  in  the  brackets  are  bounded. 

With  the  assumption  that  kh  and  hence  k'h  is  small  there  exists  >  0  s.t. 


then 


sin  kUi  =  k'h  ^1  —  -  ±  . .  >  D^k^h, 


(  - 

'■E 

>=i 


h  53  d'G.t 

i=i 


2\  5 


hi 


/  n  /i-l 

V=‘  \'=’  > 


2\  2 


and  with  the  previous  inequalities  we  obtain 


3  ^DiD 


-2 


D^hz  Z?3 


By  similar  computation  we  can  show  that  for  any  1,1  <  I  <  n 

i=i 


||//j.||  +  niax(d'//.„|  < 


D 

hVd 


.  jk'h  U  +  l)k'h 
iin  — —  cos - - - 


hence 
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where  D  =  D3  +  D^. 

Returning  now  to  eqs  (3.31)  and  (3.33), 


h^k 
<  kCD 

From  the  Taylor  scries  expansion  (3.21 )  we  see  llial 


<  f  mp  |d'//.„l  +  l|7/j)|r,l, 

\l</<n  / 

(^)  -«i- 


k‘  k^h^  ik'^h^ 

I  ”  ^  6  ~  640  ■ 

is  bounded  for  sufficiently  small  kh.  Hence  there  exists  a  constant  E  not  depending  on  h 
and  k  s.t. 

kli  <  (3.38) 

We  then  have 

|v|i  =  |u  +  c|i  <  (1  +  •Efc)|«|i, 
hence  there  exists,  for  sufficiently  large  k,  a  constant  F  s.t. 

I«li  > 

and  left  inequality  of  the  statement  follows  from  the  definition  of  z  and  the  Gardings-type 
inequality  (2.1-1). 


To  prove  the  right  inequality  we  construct,  in  analogy  to  section  2,  a  function  Zg  for 
which  continuity  holds  with  Ck~K 
Consider  the  function 

Z(x)  =  ip{x)xvix) 

where  v’(.t)  €  C‘^(0, 1)  and 


is  a  fundamental  solution  of  the  discrete  system  eq  (3.13).  Let  z^lx)  €  V*  be  the  piecewise 
linear  interpolant  of  Z(x)on  X/,.  Again  we  ajssume  that  (p  does  not  depend  on  the  paramter 
k  and  is  selected  such  that 

S?(0)  =  S?(1)  =  V’'(1)  =  0 

and  there  exists  o  >  0  s.t. 

koli  >  a 

independently  on  k.  Then 
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Turn  to  the  estimation  of  (we  omit  the  .subscript  o  from  now  on): 

Biz,v)  =  [^'v'-k^pzv 

Jo  Jo 


=  /i  d^zdh’  -  -rhY^  (Zj.i  +  4zj  +  )  v. 


j=i 


j=i 


(let  formally  :=  2„_i ). 

Summation  by  parts  then  yields 

"  /  /.2  \  1 

B(z,v)  =  -hY2^D^id^z)+  —  {Zj.i  +  4Zj  +  Zj+l)j  Vy  +  J^iZn-lVn  -  Zo%) 

The  term  outside  the  sum  is  0(h).  liulood,  r.,  =  0  and 

/,2 

=  Si>(  1 )  -  hv"'(  1 )  +  -J9"(  1 )  +  0{h^). 

Consequently,  since  I )  =  1 )  =  0,  we  have  =  />“  Vn-itt'n-i  =  0(h).  Hence, 

omitting  the  terms  0(h), 

B(z,  v)  =  -/i  ^  I^D^d^z)  +  Y  ('>-1  +  ">• 

For  arbitrarily  fi.xed  j  we  write  the  second  differences  as 

D^{d^z)  =  D^(d^(ifw))  =  £?•'  {^d^q})wj-\  +  ipjd^w^ 

=  D^(d^(f)Wj..\  +  2D^ifd^w  +  ipjD^(d^w) 


and  the  weighted  s\nn  as 


Zj-i  +  4Zj  +  ^j+i  =  +  4(s?te)j  +  (v7a')j-i 

=  u^J_^(q:j-h^p'J+0(h'‘))  +  4^Ujifj  +  ■WJ+^(ipJ■kh^p'j  +  0(h'^)) 

—  rj  (“  j-1  +  +  '"j+i )  +  2h^ip'jWj  +  0(h^)). 


Then,  neglecting  ell  term.>  that  are  0(h)  we  can  write 


/.■* 

D-’(d'z)  +  —  (:j-\  +  I'j  -f  *;+i )  = 
(» 


1.2 


D^d^w)  +  —  +  dJCj  +  M.'j+l) 


+  D^(d^q})uij-  \  +  2D^tpd^w. 


Since  w  has  been  selected  as  a  fundamental  solution  of  the  discrete  system,  the  expression 
in  square  brackets  vanislies. 

We  now  define  the  piecewise  linear  function  u  as  the  linear  interpolant  of  the  meshfunction 
Uk  defined  by 
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Then,  on  the  one  liaiicl. 

and  on  the  other  hand 


u(  1 )r( 1 

)-  /  u{x)v'{.r)(l.r 

1 

Jo 

<  (|u(l)|  +  ||n||)|e|i 


w  = 


n  /  1  —  1 


/)  b  (^D-' ((Pip)Wj^]  +  '2D^<p(Pu'"j 


‘=1  \  J=1 

1-1 


/  n  , 

\  «=1  \  J=1 

Making  use  of  the  smoothness  of  the  function  we  have  for  all  j 

£»J(r/V)  =  >p'yb)  +  Oih^) 

D^-\-  =  >?'((;  - l)/i)  +  0(/i) 


1 

2 


and  we  obtain 


Hull  <  b  (/ii|ir|  (||v"|U  +  (lly*'l|.x.  +  2||^'|k.  +  0(h)))^f 
1=1 

where  the  function  lo  =  A'“'  sin  k'x  can  be  estimated  by 


and  the  term  0(h)  does  not  depend  on 

By  similar  estimates  for  |w(l)|  we  conclude  that  for  sufficiently  small  b  there  exists  a 
constant  C  with 

(||w||  +  |t<(1)|)  <  p 

It  then  follows  that 

Veer,:  |r?(-,r)|<  ||t;|j 

and  the  proof  is  completed. 

Remark  11:  We  rccapif iilale  that,  for  /  e  1).  both  approximability  (theorem  2) 
and  the  discrete  stability  condition  hold  iiiuler  tin'  assumption  the  hk  is  sufficiently  small. 
It  then  follows  from  a  fnndamental  theorem  [B.A.  p.lKT]  that  the  FE-solution  exists  and 
is  uniquely  determined.  emphasi/.e  that  the  latter  stability  result  is  thus  obtained  by 
restricting  the  magnitude  of  /lA'  only  (compare  to  the  more  severe  restriction  of  hk^  in 
Theorem  'll). 

3.5  A  preasymptotic  error  estimate 

In  this  subsection,  an  error  estiiriate  will  be  given  that  is  suited  to  botmci  the  error  at  finite 
range  also  for  high  wavenumbers  k.  First  we  have  to  prove: 
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Lemma  3  Let  ujt  €  V/,  b<  the  finite  element  solution  to  the  variational  problem  (2.10) 
for  given  data  f  €  L^(0, 1 ). 

Then,  if  h  is  small  s.t.  Uk  <  1  ,  there  exists  a  constant  C  not  depending  on  h  and  k  s.t. 

Proof:  Since  is  piecewise  linear,  we  liave 

li'4ii  =  ('‘EC'''"/')’)’- 

Write  xth  :=  in  terms  of  the  discrete  Green’s  function  as 

n 

Ui  =  h  y  ~  GijVj 
J=» 


then 

n 

rf'w  =  hY,d'<^jrj 
j=i 

and 

lc/'«l  <  ||d’G||||r|l.  (3.39) 

with 

IHI  = 

The  mesh  function  r/,  is  related  to  the  function  /  £  L^(0, 1)  by  eq  (3.15)  from  which  it  is 
easy  to  see  thal  there  exists  a  constant  Ci  s.t. 

ikii  <  c,h^m- 


The  derivatives  of  the  Green’s  function  are  -  cf.  eqs  (3.3C,  3.37)  - 

cos  (^^{2i  —  1  ( >1  sin  k'si  +  cos  k'si)  i  <  I 

sin  k'si  ^.4  cos  (^(2?  -  1))  -  sin  ^^(2/  -  l)j^  i  >  I 

Obiously  h^  |d'6'.;|  is  bounded  provided  that  lik'  <  o  <  n/2.  From  the  Taylor  series 
expansion  of  hk'.  eq  (3.21).  we  conclude  that  .such  n  exists  for  sufficiently  (say,  hk  <  1) 
small  Ilk. 

Hence  there  is  a  constant  ('2  s.t. 


(TG.i  = 


cos  ^ 


V;'.j  ; 
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Then  also  ,, 

V. ;  IK'GII  <  ^ 

and  the  statement  follows  from  eq  (3.39)  with  C  =  CjCj.  The  proof  is  completed. 

We  proceed  to  the  proof  of  the  error  estimate: 

Theorem  5  Let  ti  6  //^(0,1)  be  the  exact  solution  of  the  variational  problem  (2.10)  with 
data  f  €  X^(0, 1)  and  let  uje  ^  5>,[0, 1(  be  the  finite  element  solution  of  (2.10). 

Then  if  the  stepsize  h  is  such  that  hk  <  1  the  error  estimate 

holds  with  constant  C  not  dejxnding  on  h  and  k. 

Proof;  Let  uj  e  Vh  =  ■StJO,  1(  be  the  interpolant  of  u  and  define  r  €  V),  by 


2  :=  it/f  -  «/. 


From 


V?'  6  \'h  '■  B(ii.i')  =  8(ujf,v) 
and  since  B  is  sosquiliiicar  wc  have 

B(u  -  iii,v)  =  B(z,  v). 

On  the  other  hand,  for  v  €  V/,: 

Vi :  f  (v  -  Uj )'  v'  =  {(n  -  «/)  -  f  («-«/)  u"  =  0 

since  («  -  w/jl.v^  =  0  and  =  0  and  therefore 

B(u- uj,v)  =  f  [u- uj)v. 

Jo 

Hence  z  is  a  solution  of 

Vv  €  Vh  :  B(z,  t>)  =  k^  (u  -  ui,v) 
and  from  Lemma  3  we  have  the  estimate 

||2'|l<CA-2||«-n/||. 

Then 

|c|i  =  |«  -  M/<|t  =  |«  -  W/  +  H/  -  U/eh 

<  |«-w/|i  +  l-li 


<  |h  -  w/li  +  -  ?tf|l 

We  now  invoke  the  approxin\ation  properties  of  the  space  Vh  from  Lemma  2  to  obtain 

M,  < 

The  statement  of  the  theorem  now  follows  from  Lemma  1. 
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3.6  Comments 

In  this  section  we  have  given  (lifTerent  proofs  of  e.\i.stence-uniqueness  for  the  FEi-solution. 
The  main  results  are  tliat 

•  the  discrete  problem  is  stalrle  provided  that  proper  restrictions  are  made  for  the 
magnitudes  of  hk  only  (theorem  4)  and 

•  the  error  of  the  finite  element  solution  can  be  controlled  by  restricting  the  magnitudes 
of  hk  and  /ifcs  (theorem  .5). 

It  had  been  shown  that  the  error  bounds  of  subsection  1  tend  towards  0  as  A;  is  increased. 
This  is  not  the  case  for  the  estimates  of  thtorem  5  since  only  boundedness  of  hk  is  assumed. 
There  is,  however,  a  close  relation  between  both  estimates.  Namely,  the  corollary  1  from 
theorem  3: 

l«  -  «/e|i  <  C\u  -  n/l,  <Ch{l  +  k)  ll/ll 

follows  also  from  theorem  5  if  the  magnitude  of  k^h  is  bounded.  In  other  words,  both  error 
estimates  lead  to  the  same  conclusion  that  the  stability  constant  Cj  does  not  depend  on 
k  if  k'^h  is  bounded. 

We  will  show  by  numerical  experiment  that  this  condition  is  also  necessary,  i.e.  the 
constant  C,  grows  w-ith  k  if  k^h  is  not  restricted. 

The  assumption  of  uniform  mesh  is  due  to  technical  necessities  of  the  proofs  for  theo¬ 
rems  4  and  5.  All  statements  of  this  .section  should  hold  for  nonuniform  mesh  as  well. 

4  Numerical  Evaluation 

The  first  and  obvious  purpo.se  of  the  numerical  evaluation  lies  in  the  illustration  and 
application  of  the  theoretical  results  by  computational  experiment.  Beyond  this,  we  will 
draw  a  qualitative  conclusion  concerning  the  a,ssumptions  of  some  propositions  of  the 
previous  section. 

Throughout  this  section,  we  will  present  FE-solutions  to  the  variational  problem  (2.10) 
with  constant  right  hand  side  J{x)  =  —1  on  uniform  mesh. 

4.1  Error  of  the  best  approximation 

Consider  in  Fig.  1  the  errors  e„  of  the  best  approximation  (interpolant)  computed  for 
different  k  and  h  such  that  0.2  <  hk  <  2,  plotted  in  log-log-scale. 

As  predicted  by  theorem  2,  all  error  curves  decrease  with  constant  slope  of  —1  in  the 
log-log-plot  (the  theoretical  rate  of  convergence  being  0(h)). 

The  inequality  of  the  theorem  gives,  however,  a  crude  upper  bound  for  the  relative 
error 


in  the  ca,so  that  /.•  is  large  and  |(//||  is  nol  bounded  from  below  independently  on  k.  For 
w  €  //^(0. 1),  an  estimale  is  obtained  from  Lemma  2  as 
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In  our  example  the  relation 

ll«'ll  V3-2cosA-  +  :i;i^i2]p^j 

is  easily  computed,  hence  for  sufficiently  large  k  we  can  by 


predict  the  meshsizc  needed  for  approximation  with  a  given  tolerance  £  =  ia- 
Let,  for  example,  f  =  0.1  be  a  maximal  tolerance,  then 


(4.2) 


is  the  ’’rule  of  the  ihuinb*’  for  the  number  of  elements.  As  table  1  shows,  this  rule  works 
well  for  large  k. 


Table  1:  Numlter  of  elements  needed  for  a  relative  error  of  interpolation  less  than  0.1: 
number  obtained  from  numerical  experintent  compared  to  bound  computed  from  eq  (4.2), 


k 

2  10  40  100 

n,  computed  from  eq  (4.2) 
Of  measured  from  Fig-  1 

Consider  now  the  results  plotted  in  Fig.  2.  Clearly  the  relative  error  of  interpolation 
cannot  exceed  100%.  From  the  plots  we  observe  that  for  each  wavenumber  k  the  error 
stays  at  100%  on  coarse  mesh  and  starts  to  decrease  at  a  certain  meshsize.  We  are  inter¬ 
ested  in  the  point  where  the  descend  starts.  More  precisely,  we  seek  the  critical  number 
of  degrees  of  freedom  according  to  the  following  definition: 

Define  -  for  any  fixed  k  and  /  -  the  critical  iiimibcr  of  degrees  of  freedom  (DOF)  No{k) 
as  the  minimal  number  N(k.f)  of  DOF  for  which 

1.  c(v.k]  <  1  and 

2.  e(ii,k)  is  monotone  decreasing  w.r.  to  ti 
for  all  n  >  N(k.f). 


For  the  best  approximation,  the  critical  number  of  DOF  is  determined  by  the  rule  that 
the  stepwidth  of  interpolation  by  piecewise  linears  should  be  smaller  than  one  half  of  the 
wavelength  of  the  exact  solution: 

Ilk  <  TT. 

In  Fig.  3,  the  critical  point  n*  computed  from 


no 


A' 

.5!-. 


(4.3) 
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is  plotted  for  difforont  k.  Tlio  iMcdicted  critiral  iiinnlier  of  DOF  is  close  to  the  actual  one 
in  all  cases. 

Finally  wc  wliisli  to  see  e.xperinieiitally  llmt  the  error  of  the  best  approximation  is 
controlled  by  bounding  the  magailnde  lik.  To  this  end,  consider  the  fat  line  plotted 
in  Fig.  3.  In  the  error  curves  for  the  different  k,  this  line  connects  the  points  that  are 
computed  from  hk  =  const.  =  0.2  .  As  we  see.  this  line  does  neither  increase  nor  decrease 
significantly  with  the  change  of  k.  For  more  detailed  ob.servation,  the  relative  error  of  the 
best  approximation,  computed  for  all  integer  k  from  1  to  500  and  for  hk  =  0.1,  is  plotted 
in  Fig.  4.  The  error  oscillates  with  decaying  amplitude  around  the  horizontal  line 

Ical,  =  0.02887. 


The  upper  estimate  from  eq  (4.1)  is 

IcJ,  <  —  =  0.03183. 

It 

The  figure  can  be  further  analyzed  as  follows:  we  find  for  the  relative  error  the  expression 
(t  =  hk)-. 

if  2  sin  2k  -  4  sin  k 

|f„|,  =  — 

under  the  assumption  that  t^  and  higher  terms  of  /  can  be  neglected.  For  the  case  t  =  0.1, 
plotted  in  Fig.  3,  this  expansion  predicts  for  high  k  the  value 

lt„|,  =  ^  =  0.02886751. 

Remark  12:  In  the  one-dimensional  ca.se  one  can  by  means  of  a  Galerkin  least  squares 
method  ([Hill])  obtain  a  modified  finite  element  solution  that  is  identical  with  the  inter- 
polant  of  the  exact  solution.  Therefore  the  conclusions  drawn  above  for  the  minimal  error 
in  /f^-seminorm  hold  for  this  solution  as  well. 

4.2  Error  of  the  finite  element  solution 

Discrete  wavenumber;  Unlike  the  best  approximation,  the  FE-solution  is,  in  general, 
not  in  phase  with  the  exact  .solution.  On  uniform  mesh  this  numerical  effect  is  highlighted 
by  the  notation  of  a  "discrete  wave  number”  k'  that  governs  the  periodicity  of  the  finite 
solution.  In  other  words,  we  observe  a  phase  lag  ([llHl,  p.71],  cf.  Fig.  4)  between  the 
exact  solution  and  it's  best  approximation  on  the  one  and  the  FE-solution  on  the  other 
hand . 

The  determining  equation  for  th  '  d:-t  relo  solution  on  uniform  mesh  had  been  found 
in  subsection  3.2.  as 


where  i  =  hk  and  the  r.li.s.  is  a  rational  function  of  t. 

In  Fig.  -5  the  functions  ij\  =  -S(t)l R(t).  yj  =  cos/  and  |.(/i|  =  1  are  plotted.  We 
observe  that: 


•  at  /o  =  \/T2  the  function  i/i  reaches  absohite  value  1;  the  numerical  solution  switches 
from  the  propagating  case  to  the  decaying  ca.se; 
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•  for  fixed  fc,  the  convergence  k'  —  k  i.s  vizualized  by  cos k'h  cost  =  cos kh  as 
h  —  0.  The  curves  begin  to  deviate  significantly  at  about  hk  =  I. 


Rate  of  convergence:  In  Fig.  6  the  relative  errors  of  the  FE-solutions  for  different 
k  are  plotted.  The  meshes  are  sucli  that  the  magnitude  of  kh  is  in  the  same  range  as  in 
the  error  plot  for  the  best  approximation  in  Fig.  1.  We  observe  the  following: 

1.  The  relative  error  of  the  FE-solution  exceeds  (for  higher  k  on  relatively  coarse  mesh) 

100%. 

2.  For  low  k  (represented  by  A-  =  3  in  the  figure)  the  rate  of  convergence  is  nearly 
constant  throughout  the  region  considered,  i.e.  the  fe-solution  behaves  essentially 
like  the  best  approximation. 

3.  For  high  k,  the  relative  error  oscillates  above  100%  before  it  starts  to  descend  after 
a  critical  value  No  of  meshpoints  has  been  reached.  The  decrease  first  occurs  with  a 
rate  greater  than  -1  in  the  log-log-scalo  but  irccomes  -1  for  small  h. 

4.  Unlike  the  error  of  the  l)cst  approxintation.  the  error  of  the  FE-solution  cannot  be 
controlled  by  hounding  the  magniinde  of  hk.  The  relative  error  clearly  grows  with 
k  on  all  lines  hk  =  const. 

The  last  observation  is  further  emphasized  by  the  results  plotted  in  Fig.  7,  together  with 
table  2.  The  ’'rule  of  the  thumb"  to  place  a  certain  number  of  elements  per  wavelength 
does  obviously  not  hold  for  high  k. 

Asymptotic  stability  and  quasioptimality:  Consider  now  in  Fig.  8  the  plots  of  the 
relative  errors  of  the  FE-solution  together  with  the  relative  errors  of  the  best  approxima¬ 
tion.  This  figure  is  well  suited  to  enhance  the  quasioptimal  stability  estimate  in  corollary 
2,  section  3.2.  To  this  end.  lines  are  plotted  connecting  h  and  k  s.t. 

hk^  =  o  =  const  (4-4) 

for  a  =  2,0  =  1  and  o  =  0.5.  The  corollary  states  that  on  these  lines,  if  o  is  sufficiently 
small,  the  ratio  of  the  errors  of  the  best  approximation  and  the  FE-solution  does  not 
depend  on  k,  i.e.  the  distances  between  both  curves  in  the  log-log-plot  do  not  grow  along 
the  lines  (4.4). 

The  statement  is  viz\ialized  in  the  plot;  even  more:  we  see  that  for  the  example  con¬ 
sidered  the  stability  constant  is  close  to  I  for  sufficiently  small  a. 

In  Fig.  9  the  stability  constant  C.,  from  corollary  2.  computed  with  the  restriction 
hk^  =  1.  is  plotted  for  all  integer  k  from  1  to  200.  Obviously,  the  constant  computed  with 
constrained  hk^  does  neither  decrea.se  nor  grow  with  increasing  k  (except  for  small  k,  then 
hk  is  the  leading  member  in  the  estimate  of  theorem  5  -  cf.  comments  to  Fig.  13). 

On  the  other  hand,  it  is  easy  to  verify  from  Fig.  that  the  error  ratio  does  depend  on  k 
on  all  lines  hk'^  =  a  with  i3  <  2.  In  particular.  C'j  is  increasing  with  k  on  the  line  defined 
by  hk  =  1  (Fig.  10)  and  hk^  =  1  (Fig.  11 ). 

Preasymptotic  stability  and  error  estimate:  We  have  thus  shown  experimentally 
that  the  assumptions  of  theorem  .3  and  corollary  2  are  indeed  inherent  to  the  problem:  for 
the  ratio  of  the  FE-solutiou  error  and  the  best  approximation  to  be  bounded  it  is  necessary 
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to  restrict  the  magnitude  of  lik^.  However,  it  is  not  nece.ssary  (tliough  sufficient)  to  bound 
this  ratio  for  the  practical  purpose  of  limiting  the  error  of  the  FE-solution  at  finite  range. 
Indeed  C,  grows  with  k  on  tlie  line  of  constant  (relative)  error  of  the  FE-solution  (Fig.  12). 

According  to  theorem  ■'j,  the  relative  error  is  bounded  at  any  range  by  the  magnitudes 
of  and  hk.  This  statement  is  vizualized  by  the  results  plotted  in  Fig.  13.  Here, 
the  relative  error  has  been  computed  for  all  integer  k  from  1  to  1000  on  meshes  with 

h  =  .  We  observe  the  following: 

1.  For  low  A:  (1  <  k  <  30)  the  relative  error  decreases  rapidly  with  k.  In  this  range, 
the  FE)-solution  is  still  clo.se  to  the  best  appro.ximation  {hk^  =  5.48  for  k  =  30)  and 
hence  the  term  hk  is  the  significant  member  in  the  estimate  (3.40). 

2.  For  large  k  (k  >  100)  the  error  is  bounded  by  c  =  0.05.  The  term  h^k^  is  leading  in 
estimate  (3.40). 

Let  us  consider  how  these  effects  might  influence  the  results  of  applied  computations. 
To  this  end,  we  write  the  estimate  of  theorem  5  in  the  form 

kl,  <  {n  +  C(l+kW)\\f\\  (4.5) 

with  the  "rule  of  the  thumb" 

hk 

—  =  o. 

In  most  practical  computations  with  low  (k  <  10)  wavenumbers,  intuitively  a  good  reso¬ 
lution  (like  a  =  0.1,  i.o.  20  elements  per  wavelength)  is  choosen.  In  this  case,  =  0.01 
and  ka"^  =  0.1;  both  terms  in  the  estimate  (4.5)  are  of  the  same  magnitude  and  hence  the 
phase  lag  does  not  affect  the  error  significantly.  In  other  words,  no  negative  effects  will 
be  observed  in  benchmark  tests.  However,  for  high  wavenumber  (say,  k  =  100)  the  second 
member  equals  1  for  the  same  resolution  o  =  0.1  and  hence  is  prevalent  in  the  estimate. 

These  effects  become  much  more  visible  if,  for  cost  reduction  of  the  computations,  there 
are  choosen  lower  resolutions  like  a  =  0.2  or  a  =  0.5  (cited  «is  ’’acceptable  resolution”  or 
"limit  of  resolution”,  respectively,  in  [Hill]).  For  k  =  10,  the  magnitudes  q  =  0.2  and 
ka^  =  0.4  are  still  of  the  same  order  for  acceptable  resolution  but  differ  considerably  for 
the  limit  of  resolution  (  o  =  0.5  and  ka^  =  2.5).  For  the  latter  resolution,  both  magnitudes 
are  roughly  of  the  same  order  up  to  k  =  -1. 

For  high  wavenumber  (k  =  100)  the  .second  member  of  the  estimate  is  clearly  domi¬ 
nating  for  both  resolutions:  we  have  o  =  0.2  vs.  ka^  =  -1  and.  for  the  limit  of  resolution, 
Q  =  0.5  and  ko^  =  25. 

Finally,  we  demonstrate  t  hat  also  the  critical  number  of  DOF  for  the  FE-solution  error 
is  governed  by  th('  magnitude  of  Consider  in  Fig.  11  the  curves  of  the  relative  error 
computed  for  different  k  from  k  =  10  to  k  =  1000  and  the  |)redicted  critical  number  of 
DOF  where  the  latter  has  been  computed  from  the  formula 


(a  physical  argument  for  this  formula  will  be  given  below).  Again,  the  predicted  critical 
number  of  DOF  is  close  to  the  actual  one. 
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4.3  Summary 

We  are  now  in  a  position  to  comment  on  the  behavior  of  both  approximation  and  FE- 
solution  error  throughout  the  whole  range  of  existence  of  a  propagating  numerical  solution 
(i.e.  for  hk  <  \/l2  in  the  case  of  piecewise  linear  approximation). 

Consider  the  case  k  =  100  in  Fig.  15.  We  have  marked  on  the  abscissa  three  signlficamt 
points  for  the  meshsize  n  =  /i~'.  By  these  points  the  range  of  degrees  of  freedom  is  divided 
into  four  regions,  namely: 

1.  1  <  n.  <  Up-. 

The  mesh  is  too  coarse  to  allow  for  neither  approximation  nor  FE-solution  of  the 
Helmholtz  equation  with  given  wavenuinber  k.  The  number 


w.,  = 


is  the  critical  number  of  DOF  (”limit  of  resolution”  [HHl])  for  approximability. 


2.  Up  <  n  <  Nq-. 

The  relative  error  of  approximation  is  smaller  than  100  %  but  the  relative  error 
of  the  FE-solution  is  still  above  this  range.  Though  we  have  approximability  and 
stability,  the  stability  constant  is  too  large  to  bound  the  error. 


3.  Np  <n<  N,: 

Both  the  FE-solution  error  and  the  approximation  error  are  in  the  range  of  con¬ 
vergence.  In  the  error  estimate  (3.^10)  the  magnitude  is  the  leading  member. 
A  considerable  phase  lag  is  present  between  the  exact  and  the  FE-solution.  The 
stability  constant 

_  |m-  M/ell 

*  inf|«-u|i 

depends  still  on  k  but  is  '’under  control”  since  the  magnitude  of  hk^^^  is  bounded. 
With  the  leading  member  of  the  estimate  being  O(h^)  (for  any  fixed  k),  the  rate 
of  convergence  of  the  FE-solution  is  higher  than  the  rate  of  convergence  of  the  best 
approximation.  The  FE-error  curve  descends  towards  the  line  of  the  optimal  error. 


4.  71  >  N,: 

The  critical  number  N,  has  been  computed  from  the  relation  hk^  =  1  (cf.  eq  (4.4) 
and  related  comments).  The  stability  constant  Cg  does  not  depend  on  h  and  k,  the 
magnitude  hk  is  leading  in  e.stimate  (3.10).  Both  the  FE-solution  error  and  the  error 
of  best  approximation  have  the  same  rate  of  convergence  0(h.),  i.e.  the  statement  of 
quasioptimality  (corollary  2)  holds. 


Concluding  this  subsection  we  give  the  argument  for  the  computation  of  the  critical 
number  of  DOF  for  the  FE-solution  eq  (4.6).  .'Vssume  that  the  solutions  are  given  by 
u  =  s'lnkx  and  «/,  =  =  sinA’'.x-/,  and  consider  the  error  in  the  £oo-norm. 

Then,  if  the  phase  lag  k  -  //  is  smaller  than  f ,  the  maximal  difference  of  amplitudes 
I  sin  kxh  —  sin  occurs  at  the  end  of  the  interval  [0, 1].  Since  ||  sin  A:a;||oo  =  1  we  require 
for  l|e||oo  <  ||sin/c.r||“’||ii  -  u/cWoc- 


k  +  k' 

.  k-k' 

cos - 

sin  — ^ — 

2 

2 

<  1- 


I  sin  k  -  sin  k'\  =  2 
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This  inequality  certainly  holds  if 


or,  equivalently, 

/,  _  /,'  <  -  I. 

"  .5 

With  this,  eq.  (4.6)  follows  from  the  Taylor  e.xpansion  eq  (3.21). 

5  Conclusions 

The  numerical  solut  ion  of  the  llclmholiz  o((uation  with  the  h- version  of  the  FEM  is  studied 
on  a  one-dimensional  model  problem.  Following  the  proof  of  new  analytical  statements, 
the  investigation  is  completed  with  the  results  of  computational  experiments. 

While  it  is  evident  from  the  oscillatory  character  of  the  exact  solution  that  the  meshsize 
h  has  to  be  adapted  to  the  magnitude  of  the  wavenumber  k  it  is  not  obvious  how  exactly 
this  adaption  should  be  properly  designed.  This  question  is  the  starting  point  and  the 
practical  motivation  of  the  present  investigation. 

On  the  one  hand,  "rules  of  the  thumb"  restricting  the  product  hk  had  reportedly  failed 
for  high  k.  On  the  other  hand,  the  restriction  of  k^h  assumed  in  existing  proofs  of  asymp¬ 
totic  stability  and  convergence  in  the  analytical  literature  are  practically  inapplicable  in 
the  very  case  of  high  wavenumbers. 

The  results  of  the  present  study  -  confined  to  the  case  of  uniform  mesh  -  reveal  that: 

•  the  finite  element  solution  is  stable  given  only  restrictions  on  the  magnitude  of  hk\ 

•  in  the  preasymplotic  range,  (he  error  of  the  finite  element  solution  is  governed  by 
the  term  li^k^  and  hence  can  be  controlled  restricting  this  magnitude; 

•  the  Babuska-Brezzi  stability  constant  is  of  order  both  in  the  infinite-dimensional 
and  the  finite-dimensional  level; 

•  the  restriction  of  hk^  is  indeed  necessary  for  quasioptimality  of  the  finite  element 
solution  w.  r.  to  k. 

In  physical  terms,  if  hk^  is  small,  then  the  FE-solution  is  in  the  asymptotic  range 
of  convergence  where  it  is  close  to  the  interpolant  of  the  exact  solution  and  hence  is 
quasioptimal,  i.e.  the  FE-error  is  proportinal  (independently  of  k)  to  the  interpolation 
error. 

In  the  preasymplotic  range,  the  difference  between  the  FE-solution  and  the  interpolant 
(the  phase  lag  of  the  FE-solution)  is  the  prevalent  part  of  the  FE-error.  To  bound  this 
error  it  is  both  necessary  and  sufficient  to  restrict  the  magnitude  li^k^. 

Referring  to  the  originally  posed  question  we  see  that  the  answer  for  the  proper  choice 
of  the  meshwidth  lies  ”in  the  middle"  (between  hk  and  hk^).  Consequently,  for  large  k 
there  still  has  to  be  chosen  a  quite  fine  mesh  (egs..  h  =  10“^  for  k  =  100)  if  the  h-version 
of  the  FEM  is  applied. 

In  part  II  of  this  paper,  results  will  be  pre.sented  for  the  h-p-version.  Following  these 
conclusions  we  have  to  investigate  what  can  be  gained  on  the  global  level  (in  terms  of  the 
number  of  DOF)  by  investing  locally  (in  terms  of  the  order  of  elemental  approximation). 


F.  Ihlenburg,  I.  Babuika,  Finite  Element  Solution  to  the  Helmholtz  Equation  ...  34 


Further  research  will  be  directed  to  the  generalization  of  the  results  presented  herein 
to  higher-dimensional  ca.ses  and  to  applied  problems  of  fluid-structure  interaction. 
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Figure  1:  Relative  error  in  //•-seminorm  for  k  =  2. 10.40  and  k  —  100. 


Figure  2:  Relative  error  of  the  best  appro.ximalion  in  //’-seminorm  and  predicted  critical 
numbers  of  DOF  for  k  =  10.  A;  =  40,  A-  =  100  and  k  =  400 
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Figure  5;  Functions  Viit)  -  <^ost 


Figure  6:  Reliitivo  ('i  ror  in  11 1-scniinorm:  Finite  clonient  solutions  for  Ic  =  3,k  =  10, 
and  k  =  100 
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Figure  9:  Stability  conslaiil  C,,  in  // '-stMiiinorni  roinputccl  with  constraint  hk^  =  1  for 
k  =  J.200. 1. 


Figure  11:  St;il)ilily  consUml  C\  in  // ‘-sc'ininorni  roni|)iil.e(l  with  constraint  =  1  for 
k  =  1,500,1. 
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Figure  12:  Rol.itivc  error  in  (I  l-semiiiorm:  Finite  element  solutions  for  k  =  10  and  k  =  100 


Figure  13:  Relative  error  of  FE-soliilion  in  //'  seminorm  computed  with  constraint 
hk^!'^  =  1  for  /.’  =  1, 1000. 1. 
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Figure  1-1:  Uolativo  error  of  llic  finite  element  solution  in  //* -semi norm  and  predicted 
critical  numbers  of  DOF  for  k  =  10.  A:  =  100.  A-  =  400  and  k  -  1000 


Figure  15:  Relative  errors  of  finite  element  solution  and  best  approximation  in  H 
seminorin  for  k  =  100  . 


The  Labmvtoiy  fw  NuiMrical  Analysis  is  an  integral  part  of  the  Institute  for  niysical 
Science  and  Tedinology  of  the  University  of  Maryland,  under  the  general  administration  of  the 
Director,  Institute  for  Physical  Science  and  Technology.  It  has  the  following  goab: 

•  To  conduct  researdi  in  the  mathematical  theory  and  computational  implementation  of 
numerical  anafysis  and  related  tqiics,  with  emphasis  on  the  numerical  treatment  of 
linear  and  nonlinear  differential  equations  and  problems  in  linear  and  nonlinear  algebra. 

•  To  help  bridge  gaps  between  computational  directions  in  engineering,  physics,  etc.,  and 
those  in  the  mathematical  community. 

•  To  provide  a  limited  consulting  service  in  all  areas  of  numerical  mathematics  to  the 
University  as  a  whole,  and  also  to  ^emment  agencies  and  industries  in  the  State  of 
Mar^nd  and  the  Washington  Metropolitan  area. 

•  To  assist  with  the  education  of  numerical  analysts,  especialfy  at  the  postdoctoral  level, 
in  conjunction  with  the  Interdisciplinary  Applied  Mathematics  Program  and  the 
programs  of  the  Mathematics  and  Computer  Sdence  Dq>artments.  This  iiududes  active 
collaboration  with  government  agencies  such  as  the  National  Institute  of  Standards  and 
Technology. 

•  To  be  an  international  center  of  study  and  research  for  foreign  students  in  numerical 
mathematics  who  are  supported  by  foreign  governments  or  exchange  agencies 
(Fulbri^t,  etc.). 

Further  information  may  be  obtained  from  Professor  I.  Babudka,Chairman,  Laboratory  for 
Numerical  Analysis,  Institute  for  Physical  Science  and  Tedmology,  University  of  Maryland,  College 
Park,  Maryland  20742-2431. 


